Posterior sampling ======================== .. index:: rsamplers.rsample.sample The sampling function ----------------------- .. include:: rsamplers_sample.rst .. index:: rsamplers.rwmh Algorithm : Random-walk Metropolis Hastings --------------------------------------------- .. include:: rsamplers_rwmh.rst .. index:: rsamplers.imh Algorithm : Independent Metropolis Hastings ------------------------------------------------ .. include:: rsamplers_imh.rst .. index:: rsamplers.apt Algorithm : Adaptive parallel tempering ------------------------------------------ .. include:: rsamplers_apt.rst .. index:: rsamplers.slice Algorithm : Slice Sampler ------------------------------ .. include:: rsamplers_slice.rst .. index:: rsamplers.usrsmplr Algorithm : usrsmplr Sampler ------------------------------ .. include:: rsamplers_usrsmplr.rst Delayed acceptance: surrogate-accelerated chains ------------------------------------------------- The Metropolis-Hastings samplers (``rwmh``, ``imh``) accept a ``da_surrogate`` option: a cheap approximation of the log posterior (e.g. ``@(x)predict(s, x)`` from a surrogate fitted by ``emulate``) that screens every proposal before the true posterior is evaluated. Only proposals that survive the screen pay for a true evaluation, and a second accept/reject step keeps the chain's stationary distribution **exactly** the true posterior. See the :doc:`Surrogates chapter <../GSAandUQ/surrogates>` for the full story and a worked fs2000 example (~70% of true evaluations skipped, ~3x wall time). Checkpointing and resuming chains ---------------------------------- Long runs should not die with a crash. The Metropolis-Hastings samplers accept three persistence options:: opts.store_file = 'mychain.mat'; % checkpoint file (MAT v7.3 = HDF5) opts.store_every = 500; % checkpoint frequency (iterations) opts.resume = false; % restart from the checkpoint With ``store_file`` set, the complete chain state -- draws so far, current position, rng state, adaptive-tuning state, and the delayed-acceptance surrogate value when active -- is written every ``store_every`` iterations (atomically: a crash mid-write cannot corrupt the previous checkpoint). Setting ``resume = true`` restarts from the checkpoint, and because the rng and tuning states are restored exactly, the continued chain is **bit-for-bit** the chain that would have run uninterrupted. A completed run can also be *extended*: resume with a larger ``N`` and the chain picks up where it stopped. With ``nchain > 1``, ``_chain`` is appended to the file name per chain. The tempered ``apt`` and ``slice`` samplers do not support checkpointing yet and error early rather than silently ignoring the option. User-defined Algorithms: Another approach ------------------------------------------- Due to the fact that the RISE toolbox is written in a modular manner, a user can apply his own sampling algorithm. However, if he wants RISE to process the draws, the output of the sampling should be organized in a way that RISE understands. Here is an example in which we write a wrapping function **wrapped\_dramrun** around the delayed-rejection adaptive Metropolis (DRAM) sampler available `here `_ . .. code-block:: matlab :linenos: function r=wrapped_dramrun(m,x0,SIG,N) % wrapped_dramrun : RISE wrapper for the DRAM code [~,lb,ub,x0_,vcov]=pull_objective(m); [model,data,param,opts]=dramrun_inputs(); % call to the DRAM code [results, chain] = dramrun(model, data, param, opts); r=results2rise(); function [model,data,param,opts]=dramrun_inputs() if isempty(SIG),SIG=vcov; end if isempty(SIG),x0=x0_; end % RISE assumes column vector, DRAM assumes row vector param = struct(); param.par0 = x0.'; param.bounds=[lb,ub].'; opts=struct(); opts.qcov = SIG; opts.nsimu=N;% Number of samples data = {}; model=struct(); model.ssfun = @(p)m.routines.likelihood(p.',m); model.priorfun = @(p)log_prior_density(m,p.'); end function r=results2rise() draws=chain.'; r=struct(); for id=1:N r.pop(id)=struct('x',draws(:,id),'f',nan); end r.stats.accept_ratio=results.accepted; r.stats.funevals=nan; r.stats.time_elapsed=nan; r.SIG=results.qcov; % covariance matrix now hidden r.c=results.adascale^2; end end As the code shows, there are three main articulations of the program : #. Obtain from RISE and from the user the inputs that DRAM requires #. Run DRAM #. Transform the output of DRAM so that it can be further processed by RISE Obtaining the DRAM inputs ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Some samplers can be run directly using the log-posterior kernel, which can be obtain in RISE through the function ``pull_objective``. The DRAM code on the other hand, requires one function to evaluate log-prior and a separate function to evaluate the log-likelihood. In RISE, we can get a handle to the function that evaluates the prior by using ``log_prior_density``. We need to remember though that RISE stores the parameters vertically while DRAM stores them horizontally. Therefore, before evaluating the prior, we need to transpose the vector of parameters. .. code-block:: matlab prior = @(theta) log_prior_density(model, theta.'); Where "model" is the ``rise_model`` object. Likewise we can obtain the function computing the likelihood as follows: .. code-block:: matlab likelihood = @(theta) - model.routines.likelihood(theta', model); The other elements required by the DRAM code are easier to collect. These include : the starting vector (x0), the number of draws (n), the covariance matrix (qcov), additional inputs to the likelihood function if any (data). Running DRAM ~~~~~~~~~~~~~~~~~~~ With the inputs collected, we can then go ahead and run the DRAM sampler. .. code-block:: matlab [results, chain, s2chain] = DRAM code(model, data, param, options); Transforming the DRAM output ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ If we want to process the output of DRAM through RISE, we need to transform the DRAM output in a form that RISE expects. RISE expects each vector to be stored in a structure alongside the value of the posterior kernel. Unfortunately DRAM does not return the value of the posterior kernel. This is why we just set to "nan" all elements that cannot be obtained from the DRAM code. Processing Posterior draws : The mcmc class ============================================= .. toctree:: mcmcHelper .. # include:: mcmcHelper.rst causes problems of "duplicate label" Visualizing posteriors (and priors) ------------------------------------ See :ref:`Visualizing priors and posteriors`. Marginal Data Density computation : the mdd class ==================================================== .. toctree:: marginalDataDensity .. # include:: marginalDataDensity.rst causes problems of "duplicate label"